Introduction

In this lab we’ll introduce the add-on package tmap, which can be used to produce much nicer maps than the default mapping functions in R. We’ll also look briefly at leaflet, a package for making interactive maps. Before starting the lab, make sure both of these packages are installed. Other useful ones include mapsf highcharter, mapview and mapdeck.

Remember to use RStudio’s script editor (top left panel) to enter your commands and then copy-paste these to the console when you are ready. This will keep a copy of the commands that you can save and return to.

Before starting, remember to create a working directory (e.g. module15). Next download the files for today’s lab from the Canvas page and move them to this directory (unzip the compressed files). You’ll need the following files:

  • Climate dataset for Western North America: WNAclimate.csv
  • A shapefile of UTA light rail routes: LightRail_UTA.shp
  • A shapefile of UTA light rail stations: LightRailStations_UTA.shp
  • A shapefile of Salt Lake City census tracts: slctracts.shp
  • A LandSat8 scene for part of California: rs.zip
  • A shapefile of California places: ca_places.shp

Note for the shapefiles, you will need to download and decompress the zip file for each one to get the shapefile and the associated extra files. If you are not sure about this, please let me know.

Now start RStudio, and change your working directory from the [Session] menu > [Set working directory] > [Choose directory…]. Remember to check that R can see the files by running the list.files() command.

Once you have done this, let’s load all the files we will need:

library(sf)
library(terra)

sf_use_s2(FALSE)

tmap

tmap makes thematic maps. It works in a very similar way to ggplot2, by building a series of layers and map geometries and elements. In general, we start by using tm_shape() to identify the spatial object to be used, and add geometries are added, including filled polygons, borders and symbols. Finally, we can add legends, scale bars, etc. We’ll look at three examples of working with this package here, but there are a lot more details and examples on the [TMap website][tmapid].

Example 1: Salt Lake County

We’ll first make a map with a mixture of vector layers, including census tracts, light rail tracks and stations. Let’s start by loading these:

## Census tracts for Salt Lake County with population density
tracts <- st_read("./slc_tract/slc_tract_2015.shp", quiet = TRUE)
## Salt Lake light rail tracks
lightrail <- st_read("./LightRail_UTA/LightRail_UTA.shp", quiet = TRUE)
## Salt Lake light rail stations
stations <- st_read("./LightRailStations_UTA/LightRailStations_UTA.shp", quiet = TRUE)

If we check the CRS, you’ll see that tracts has a different one, so will need to be reprojected:

st_crs(tracts)$epsg
## [1] 4269
st_crs(lightrail)$epsg
## [1] 26912
tracts <- st_transform(tracts, st_crs(lightrail))

Now let’s start mapping things out with tmap. First, let’s make a simple map showing the polygon outlines using tm_borders(). This is the standard tmap method: first define the object that you want to plot with tm_shape, then define what you want to plot.

library(tmap)

tm_shape(tracts) + 
  tm_borders()

The function tm_fill() will then fill these using one of the variables in the tract data set. We’ll use the population density (density). Note that this automatically adds a legend within the frame of the figure:

tm_shape(tracts) + 
  tm_fill("density") +
  tm_borders()

The color scale can be changed by setting the palette argument in tm_fill(). This includes ColorBrewer and viridis scales and there are a set of methods todefine the intervals in the color palette. For example, to use the ‘Greens’ palette with percentile breaks. We also change the color of the tract boundaries to lightgray:

tm_shape(tracts) + 
  tm_fill("density", palette = "Greens", style = "quantile") +
  tm_borders("lightgray")

You can reverse the color scale by prepending a -. This uses a reversed magma palette from the viridis set, with a continuous palette:

tm_shape(tracts) + 
  tm_fill("density", palette = "-magma", style = "cont") +
  tm_borders("lightgray")

If you want to see the full set of palettes that you can use, install the tmaptools package and run the following code (you may also need to install shinyjs):

tmaptools::palette_explorer()

Let’s now add another layer. We’ll overlay the light rail track on this map. As this is a different object in R (lightrail), we first need to use tm_shape to indicate that we are using it, then add tm_lines() to show the routes. We’ll make the lines a little thicker with lwd and change the type to dashed.

tm_shape(tracts) + 
  tm_fill("density", palette = "Greens", style = "quantile") +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 2, lty = 'dashed', col = "darkorange")

Note that if you use a varible name for the color, it will thematically map the lines (I’ve removed the polygon fill to make this clear):

tm_shape(tracts) + 
  #tm_fill("density", palette = "Greens", style = "quantile") +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 4, col = "ROUTE")

We can now add the stations (using tm_shape() to select the correct object to plot):

tm_shape(tracts) + 
  tm_fill("density", palette = "Greens", style = "quantile") +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 2) +
  tm_shape(stations) +
  tm_dots(size = 0.25, shape = 23)

Now let’s add some details to the map. We’ll add a graticule before the polygons are plotted, a compass and a scalebar. We’ll also use the tm_layout function to add a title and move the legend to a new position. (This function has a lot of options for changing your map - it’s well worth checking the help page.)

tm_shape(tracts) + 
  tm_graticules(col = "lightgray") + 
  tm_fill("density", title = "Popn density", palette = "Greens", style = "quantile") +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 2) +
  tm_shape(stations) +
  tm_dots(size = 0.25, shape = 23) +
  tm_compass(position = c("left", "bottom")) +
  tm_scale_bar(position = c("right", "top")) +
  tm_layout(main.title = "Salt Lake County Light Rail", 
            legend.outside = TRUE,
            legend.outside.position = c("left"))

The map can be clipped to smaller regions by using the bbox argument in the first tm_shape() function. In this code, we’ll extract a set of census tracts that correspond to downtown Salt Lake. We’ll then use the the bounding box of these to crop the map (I’ve removed the fill). We’ll also add the station name using tm_text():

tracts_sub = tracts %>%
  dplyr::filter(TRACTCE %in% c(102500, 102600, 114000))

tm_shape(tracts, bbox = st_bbox(tracts_sub)) + 
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 2) +
  tm_shape(stations) +
  tm_dots(size = 0.25, shape = 23) +
  tm_text("STATIONNAM", ymod = -1, bg.color = "white", size = 0.8)

Finally, tmap allows you to make interactive maps fairly simply with the tmap_mode() function. Setting this to view will make an interactive map, to plot will make a static map (like the ones we have made so far). Here we’ll remake the full map as an interactive map:

## Set interactive
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(tracts) + 
  tm_fill("density", title = "Popn density", palette = "Greens", style = "quantile", id = "TRACTCE") +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 2) +
  tm_shape(stations) +
  tm_dots(size = 0.25, shape = 23)
## Symbol shapes other than circles or icons are not supported in view mode.

We’ll reset back to static maps for the rest of this exercise:

tmap_mode("plot")
## tmap mode set to plotting

Example 2: Western North America climate

In this example, we’ll map out the climate data from western North America. We’ll go a bit faster in this example and not explain every step. First we’ll need to convert the wna_climate data frame to an sf object:

wna_climate <- read.csv("./WNAclimate.csv")
wna_climate <- st_as_sf(wna_climate, 
                        coords = c("LONDD", "LATDD"),
                        crs = 4326)

Make a quick thematic map of average july temperature:

tm_shape(wna_climate) +
  tm_symbols(col = "Jul_Tmp")

Next, we’ll download country outlines and river centerlines from [Natural Earth][natEarthID] using the rnaturalearth package (you will need to install this). This will download layers from the Natural Earth website - there are a variety of these available, including country and state boundaries, physical features (rivers, lakes, etc) and cultural features (place names, roads, etc). There are three levels of resolution: 1:10m, 1:50m and 1:110m. To download a layer, you need to specify the resolution (scale) and the name of the layer. The full set of available layers can be found [here][rneid]

library(rnaturalearth)

countries50 <- ne_download(scale = 50, type = "admin_0_countries")
rivers50 <- ne_download(scale = 50, type = "rivers_lake_centerlines", category = "physical")

Now we can put this all together as a series of layers. We’ll make the color palette continuous and use a red-blue palette.

tm_shape(countries50, bbox = st_bbox(wna_climate)) +
  tm_borders() +
  tm_shape(rivers50) +
  tm_lines("lightblue", lwd = 2) +
  tm_shape(wna_climate) +
  tm_symbols(col = "Jul_Tmp", palette = "-RdBu", alpha = 0.75,
             style = "cont", title.col = "degC")
## Warning: The shape rivers50 contains empty units.

Finally, we’ll make two maps (july temperature and annual precipitation), and use tmap_arrange() to make a two panel plot:

m1 = tm_shape(countries50, bbox = st_bbox(wna_climate)) +
  tm_borders() +
  tm_shape(rivers50) +
  tm_lines("lightblue", lwd = 2) +
  tm_shape(wna_climate) +
  tm_symbols(col = "Jul_Tmp", palette = "-RdBu", alpha = 0.75,
             style = "cont", title.col = "degC")

m2 = tm_shape(countries50, bbox = st_bbox(wna_climate)) +
  tm_borders() +
  tm_shape(rivers50) +
  tm_lines("lightblue", lwd = 2) +
  tm_shape(wna_climate) +
  tm_symbols(col = "annp", palette = "BuPu", alpha = 0.75,
             style = "cont", title.col = "mm/yr")

tmap_arrange(m1, m2)
## Warning: The shape rivers50 contains empty units.
## Warning: The shape rivers50 contains empty units.
## Warning: The shape rivers50 contains empty units.
## Warning: The shape rivers50 contains empty units.

Example 3: Raster images

In the final example, we’ll use the Landsat data with tmap.

## Landsat images for California
filenames <- paste0('./rs/LC08_044034_20170614_B', 1:11, ".tif")
landsat <- rast(filenames)
## Places for California
ca_places <- st_read("./ca_places/ca_places.shp", quiet = TRUE)

The main function to map raster data is tm_raster(). We’ll remake the NDVI layer that we made in lab 12 and plot this.

b4 <- landsat[[4]]
b5 <- landsat[[5]]

ndvi <- (b5 - b4) / (b5 + b4)

plot(ndvi, col=rev(terrain.colors(10)), main = "NDVI")

tm_shape(ndvi) +
  tm_raster()

We’ll update this a little by clamping the NDVI values at a lower limit of 0, changing the color palette and making it continuous. We’ll also add the California places added as an overlay.

tm_shape(clamp(ndvi, lower = 0)) +
  tm_raster(palette = "Greens", style = "cont", title = "NDVI") +
  tm_shape(ca_places) +
  tm_borders(lwd = 2) +
  tm_layout(main.title = "Landsat 8 (2017/06/14)", 
            legend.position = c("left", "top"), 
            legend.bg.color = "white", legend.bg.alpha = 0.7)

A second function (tm_rgb()) allows the creation of color composites. We’ll use the red, green and blue bands to make the true color composite. The maximum reflectance value in these files is 1.0, so we need to set this for scaling.

tm_shape(landsat) + 
  tm_rgb(r = 4, g = 3, b = 2, max.value = 1)
## Warning: Raster values found that are outside the range [0, 1]

When we made this figure earlier, it was a lot brighter. This is because the previous function ‘stretched’ the range of values for each layer. We can mimic that here by first using the stretch() function from terra, then repeating the plot:

landsat_stretch = stretch(landsat, minv = 0, maxv = 1, 
                          minq = 0.01, maxq = 0.99)

tm_shape(landsat_stretch) + 
  tm_rgb(r = 4, g = 3, b = 2, max.value = 1)

leaflet

The tmap mode function allows you to quickly make an interacitve map. The resulting map is based on the well-known Leaflet javascript library, which offers a much greater range of flexibility in making these maps. Fortunately, there is an R library (leaflet) that allows access to the Leaflet API allowing you to build quite complex maps from R. We’ll remake the map of Salt Lake County and the light rail routes using this now. First install the library (install.packages("leaflet")), and load it (and a helper library for working with html formats):

library(leaflet)
library(htmltools)

leaflet operates in a similar way to ggplot2 and tmap in building a map by layers. Layers are added using the pipe operator that we encountered in an earlier lab (%>%). There are three basic ingredients that are required:

  • A blank map canvas from the leaflet() function
  • A set of background tiles
  • One or more markers or polygons to display

Here’s a really simple example, where we use the default tiles from OpenStreetMap and drop a marker on the University of Utah:

m <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addMarkers(lng=-111.8421, lat=40.7649, popup="The University of Utah")
m  # Print the map

You can also try changing the Provider - the source of the background map. For example, this will use ESRI map (the map is not shown here but will appear in RStudio when you run this):

m <- leaflet() %>%
  addProviderTiles("Esri.WorldStreetMap") %>%  
  addMarkers(lng=-111.8421, lat=40.7649, popup="The University of Utah")
m

ESRI imagery:

m <- leaflet() %>%
  addProviderTiles("Esri.WorldImagery") %>%  # Add default OpenStreetMap map tiles
  addMarkers(lng=-111.8421, lat=40.7649, popup="The University of Utah")
m

You can find a full set of these providers here: https://leaflet-extras.github.io/leaflet-providers/preview/

We’ll now add the census tracts and light rail sf objects. First, these need to be reprojected to WGS84 long-lat (EPSG code 4326) to work with Leaflet.

tracts_ll <- st_transform(tracts, crs = 4326)
stations_ll <- st_transform(stations, crs = 4326)
lightrail_ll <- st_transform(lightrail, crs = 4326)

Now we can add these to a dark basemap. First the tracts with addPolygons, then the routes (addPolylines) and stations (addMarkers). Note that the stations include a label taken from the station name column in the sf object. This wil appear when you hover your cursor over a marker.

leaflet() %>%
  # add a dark basemap
  addProviderTiles("CartoDB.DarkMatter") %>%
  # add the polygons of the clusters
  addPolygons(
    data = tracts_ll,
    color = "#E2E2E2",
    opacity = 1, # set the opacity of the outline
    weight = 1, # set the stroke width in pixels
    fillOpacity = 0.2 # set the fill opacity
  ) %>%
  addPolylines(
    data = lightrail_ll
  ) %>%
  addMarkers(
    data = stations_ll,
    label = ~htmlEscape(STATIONNAM)
  )

Finally, we’ll add some extra details. First, we’ll set colors for each station by the light rail line it lies on. This is a little complicated - below we:

  • Get the list of unique station names
  • Make up a color palette from the RColorBrewer package (you may need to install this)
  • Convert the station names to an integer code (so the first name alphabetically will be 1, etc)
  • Finally use the integer code to assign the matching color back to the sf object
unique(stations_ll$LINENAME)
## [1] "University"            "N/S TRAX"              "Med Center"           
## [4] "Intermodal Hub"        "Airport"               "Draper"               
## [7] "Mid-Jordan"            "West Valley"           "Sugar House Streetcar"
line_pal <- RColorBrewer::brewer.pal(9, "Set1")
line_no <- as.numeric(as.factor(stations_ll$LINENAME))

stations_ll$LINECOL <- line_pal[line_no]

Next we’ll make a popup window for each station. This will appear when the marker is clicked (as opposed to the label which appears when a cursor hovers over the marker). This has to be formatted using html. Here, we use R’s paste() function to stick together html tags with the station name, address, etc. The <b> tags make the labels bold, and the <br> tags add a line break. You can do a lot more here, for example, including images or URLs in each popup.

station_popup = paste0(
  "<b>Station: </b>",
  stations_ll$STATIONNAM,
  "<br>",
  "<b>Line Name: </b>",
  stations_ll$LINENAME,
  "<br>",
  "<b>Park n Ride: </b>",
  stations_ll$PARKNRIDE,
  "<br>",
  "<b>Address: </b>",
  stations_ll$ADDRESS
  
)

With all this in place, we can rebuild our map. The main changes here are that we use a circle marker as this has an option to change color, and include the popup information in the marker.

leaflet() %>%
  # add a dark basemap
  addProviderTiles("Esri.WorldStreetMap") %>%
  # add the polygons of the clusters
  addPolygons(
    data = tracts_ll,
    color = "#E2E2E2",
    # set the opacity of the outline
    opacity = 1,
    # set the stroke width in pixels
    weight = 1,
    # set the fill opacity
    fillOpacity = 0.2
  ) %>%
  addPolylines(
    data = lightrail_ll
  ) %>%
  addCircleMarkers(
    data = stations_ll,
    color = ~LINECOL,
    label = ~htmlEscape(STATIONNAM),
    popup = station_popup
  )

Exercise

There is no additional exercise for today, but you need to submit your R script or (preferably) a Markdown document with the code examples from the lab